#include "nnec/nnec_kernels.h"
#include "nnec/nnec_constants.h"

using namespace std;

namespace NNEC
{
    static double Integrand_G_R(double x, void * params);
    static double Integrand_G_I(double x, void * params);
    static double Integrand_G1_R(double x, void * params);
    static double Integrand_G1_I(double x, void * params);
    static double Integrand_G2_R(double x, void * params);
    static double Integrand_G2_I(double x, void * params);
    static double Integrand_G3_R(double x, void * params);
    static double Integrand_G3_I(double x, void * params);
    static double Integrand_G4_R(double x, void * params);
    static double Integrand_G4_I(double x, void * params);
    static double Integrand_Gt_R(double x, void * params);
    static double Integrand_Gt_I(double x, void * params);

    complex<double> G(double rho, double z, double z1, double z2)
    {
        double k = M_2PI;

        double g_r, g_i;
        int neval;

        INTEGRAND integrand_g;

        Integrand_Parameter p;
        p.rho = rho;
        p.z = z;
        p.k = k;

        integrand_g.parameter = &p;

        integrand_g.function = &Integrand_G_R;
        integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_r, &neval);

        integrand_g.function = &Integrand_G_I;
        integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_i, &neval);

        complex<double> value = complex<double>(g_r, g_i);

        return value;
    }

    complex<double> G1(double rho, double a, double z, double z0)
    {
        double g1_r, g1_i;
        int neval;

        double k = M_2PI;

        INTEGRAND integrand_g1;

        Integrand_Parameter p;
        p.rho = rho;
        p.a = a;
        p.k = k;
        p.z = z - z0;

        integrand_g1.parameter = &p;

        integrand_g1.function = &Integrand_G1_R;
        integration_romberg(&integrand_g1, 0.0, M_2PI, 20, 1e-10, 1e-7, &g1_r, &neval);

        integrand_g1.function = &Integrand_G1_I;
        integration_romberg(&integrand_g1, 0.0, M_2PI, 20, 1e-10, 1e-7, &g1_i, &neval);

        complex<double> value = 0.5 * complex<double>(g1_r, g1_i) / M_PI;

        return value;
    }

    complex<double> G2(double rho, double a, double z, double z0)
    {
        double g2_r, g2_i;
        int neval;

        double k = M_2PI;

        INTEGRAND integrand_g2;

        Integrand_Parameter p;
        p.rho = rho;
        p.a = a;
        p.k = k;
        p.z = z - z0;

        integrand_g2.parameter = &p;

        integrand_g2.function = &Integrand_G2_R;
        integration_romberg(&integrand_g2, 0.0, M_2PI, 20, 1e-10, 1e-7, &g2_r, &neval);

        integrand_g2.function = &Integrand_G2_I;
        integration_romberg(&integrand_g2, 0.0, M_2PI, 20, 1e-10, 1e-7, &g2_i, &neval);

        complex<double> value = 0.5 * complex<double>(g2_r, g2_i) / M_PI;

        return value;
    }

    complex<double> G3(double rho, double a, double z, double z0)
    {
        double g3_r, g3_i;
        int neval;

        double k = M_2PI;

        INTEGRAND integrand_g3;

        Integrand_Parameter p;
        p.rho = rho;
        p.a = a;
        p.k = k;
        p.z = z - z0;

        integrand_g3.parameter = &p;

        integrand_g3.function = &Integrand_G3_R;
        integration_romberg(&integrand_g3, 0.0, M_2PI, 20, 1e-10, 1e-7, &g3_r, &neval);

        integrand_g3.function = &Integrand_G3_I;
        integration_romberg(&integrand_g3, 0.0, M_2PI, 20, 1e-10, 1e-7, &g3_i, &neval);

        complex<double> value = 0.5 * complex<double>(g3_r, g3_i) / M_PI;

        return value;
    }

    complex<double> G4(double rho, double a, double z, double z1, double z2)
    {
        double k = M_2PI;

        double g_r, g_i;
        int neval;

        INTEGRAND integrand_g;

        Integrand_Parameter p;
        p.rho = rho;
        p.a = a;
        p.z = z;
        p.k = k;

        integrand_g.parameter = &p;

        integrand_g.function = &Integrand_G4_R;
        integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_r, &neval);

        integrand_g.function = & Integrand_G4_I;
        integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_i, &neval);

        complex<double> value = 0.5 * complex<double>(g_r, g_i) / M_PI;

        return value;
    }

    complex<double> G1A(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        complex<double> value;

        if(rho > a) {
            double Rp2 = rho * rho + (z - z0) * (z - z0);
            double Rp = sqrt(Rp2);
            double Rp5 = Rp * Rp2 * Rp2;
            complex<double> emjkRp = complex<double>(cos(k * Rp), -sin(k * Rp));

            value = rho * rho * complex<double>(3.0 - k * k * Rp2, 3.0 * k * Rp)
                    - 2.0 * Rp2 * complex<double>(1.0, k * Rp);
            value = emjkRp / Rp + 0.25 * a * a * emjkRp * value / Rp5;
        }
        else {
            double Ra2 = a * a + (z - z0) * (z - z0);
            double Ra = sqrt(Ra2);
            double Ra5 = Ra * Ra2 * Ra2;
            complex<double> emjkRa = complex<double>(cos(k * Ra), -sin(k * Ra));

            value = a * a * complex<double>(3.0 - k * k * Ra2, 3.0 * k * Ra)
                    - 2.0 * Ra2 * complex<double>(1.0, k * Ra);
            value = emjkRa / Ra + 0.25 * rho * rho * emjkRa * value / Ra5;
        }

        return value;
    }

    complex<double> G2A(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        complex<double> value;

        if(rho > a) {
            double Rp2 = rho * rho + (z - z0) * (z - z0);
            double Rp = sqrt(Rp2);
            double Rp5 = Rp * Rp2 * Rp2;
            complex<double> emjkRp = complex<double>(cos(k * Rp), -sin(k * Rp));

            value = emjkRp / rho / Rp + 0.125 * a * a * rho * (a * a / rho / rho + 2.0)
                    * complex<double>(3.0 - k * k * Rp2, 3.0 * k * Rp) * emjkRp / Rp5;
        }
        else if (rho < a) {
            double Ra2 = a * a + (z - z0) * (z - z0);
            double Ra = sqrt(Ra2);
            double Ra3 = Ra * Ra2;
            double Ra5 = Ra2 * Ra3;
            complex<double> emjkRa = complex<double>(cos(k * Ra), -sin(k * Ra));

            value = -0.5 * rho * complex<double>(1.0, k * Ra) * emjkRa / Ra3;
            value = value - 0.125 * rho * rho * rho
                    * complex<double>(3.0 - k * k * Ra2, 3.0 * k * Ra) * emjkRa / Ra5;
        }
        else {
            double Ra2 = a * a + (z - z0) * (z - z0);
            double Ra = sqrt(Ra2);
            double Ra5 = Ra * Ra2 * Ra2;
            complex<double> emjkRa = complex<double>(cos(k * Ra), -sin(k * Ra));

            value = a * a * complex<double>(3.0 - k * k * Ra2, 3.0 * k * Ra)
                    - 2.0 * Ra2 * complex<double>(1.0, k * Ra);
            value = 0.5 * emjkRa / a / Ra + 0.125 * a * emjkRa * value / Ra5;
        }

        return value;
    }

    complex<double> G3A(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        complex<double> value;

        if(rho > a) {
            double Rp2 = rho * rho + (z - z0) * (z - z0);
            double Rp = sqrt(Rp2);
            double Rp3 = Rp * Rp2;
            complex<double> emjkRp = complex<double>(cos(k * Rp), -sin(k * Rp));

            value = emjkRp / rho + complex<double>(0.0, 0.125 * k * a * a * rho)
                    * (a * a / rho / rho + 2.0) * complex<double>(1.0, k * Rp)
                    * emjkRp / Rp3;
        }
        else if (rho < a) {
            double Ra2 = a * a + (z - z0) * (z - z0);
            double Ra = sqrt(Ra2);
            double Ra3 = Ra * Ra2;
            complex<double> emjkRa = complex<double>(cos(k * Ra), -sin(k * Ra));

            value = complex<double>(0.0, -0.5 * k * rho) * emjkRa / Ra
                    - complex<double>(0.0, 0.125 * k * rho * rho * rho) * complex<double>(1.0, k * Ra)
                    * emjkRa / Ra3;
        }
        else {
            double Ra2 = a * a + (z - z0) * (z - z0);
            double Ra = sqrt(Ra2);
            double Ra3 = Ra * Ra2;
            complex<double> emjkRa = complex<double>(cos(k * Ra), -sin(k * Ra));

            value = 0.5 * emjkRa / a
                    - complex<double>(0.0, 0.25 * k * a) * emjkRa / Ra
                    + complex<double>(0.0, 0.125 * k * a * a * a) * complex<double>(1.0, k * Ra)
                    * emjkRa / Ra3;
        }

        return value;
    }

    complex<double> G4A(double rho, double a, double z, double z1, double z2)
    {
        double k = M_2PI;
        complex<double> c1, c2, value;
        double g_r, g_i;
        int neval;

        INTEGRAND integrand_g;
        Integrand_Parameter p;

        if(rho > a) {
            double Rp2 = rho * rho + (z - z1) * (z - z1);
            double Rp = sqrt(Rp2);
            double Rp3 = Rp * Rp2;
            c1 = 0.25 * a * a * (z - z1) * complex<double>(1.0, k * Rp)
                 * complex<double>(cos(k * Rp), -sin(k * Rp)) / Rp3;

            Rp2 = rho * rho + (z - z2) * (z - z2);
            Rp = sqrt(Rp2);
            Rp3 = Rp * Rp2;
            c2 = 0.25 * a * a * (z - z2) * complex<double>(1.0, k * Rp)
                 * complex<double>(cos(k * Rp), -sin(k * Rp)) / Rp3;

            p.rho = rho;
            p.z = z;
            p.k = k;

            integrand_g.parameter = &p;

            integrand_g.function = &Integrand_G_R;
            integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_r, &neval);

            integrand_g.function = & Integrand_G_I;
            integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_i, &neval);

            value = (1.0 - 0.25 * k * k * a * a) * complex<double>(g_r, g_i) - (c2 - c1);
        }
        else {
            double Ra2 = a * a + (z - z1) * (z - z1);
            double Ra = sqrt(Ra2);
            double Ra3 = Ra * Ra2;
            c1 = 0.25 * rho * rho * (z - z1) * complex<double>(1.0, k * Ra)
                 * complex<double>(cos(k * Ra), -sin(k * Ra)) / Ra3;

            Ra2 = a * a + (z - z2) * (z - z2);
            Ra = sqrt(Ra2);
            Ra3 = Ra * Ra2;
            c2 = 0.25 * rho * rho * (z - z2) * complex<double>(1.0, k * Ra)
                 * complex<double>(cos(k * Ra), -sin(k * Ra)) / Ra3;

            p.rho = a;
            p.z = z;
            p.k = k;

            integrand_g.parameter = &p;

            integrand_g.function = &Integrand_G_R;
            integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_r, &neval);

            integrand_g.function = & Integrand_G_I;
            integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_i, &neval);

            value = (1.0 - 0.25 * k * k * rho * rho) * complex<double>(g_r, g_i) - (c2 - c1);
        }

        return value;
    }

    complex<double> G1a(double rho, double z, double z0)
    {
        double k = M_2PI;
        double Rp = sqrt(rho * rho + (z - z0) * (z - z0));
        complex<double> value = complex<double>(cos(k * Rp), -sin(k * Rp)) / Rp;
        return value;
    }

    complex<double> G2a(double rho, double z, double z0)
    {
        double k = M_2PI;
        double Rp = sqrt(rho * rho + (z - z0) * (z - z0));
        complex<double> value = complex<double>(cos(k * Rp), -sin(k * Rp)) / rho / Rp;
        return value;
    }

    complex<double> G3a(double rho, double z, double z0)
    {
        double k = M_2PI;
        double Rp = sqrt(rho * rho + (z - z0) * (z - z0));
        complex<double> value = complex<double>(cos(k * Rp), -sin(k * Rp)) / rho;
        return value;
    }

    complex<double> G4a(double rho, double z, double z1, double z2)
    {
        double g_r, g_i;
        int neval;

        double k = M_2PI;

        INTEGRAND integrand_g;

        Integrand_Parameter p;
        p.rho = rho;
        p.k = k;
        p.z = z;

        integrand_g.parameter = &p;

        integrand_g.function = &Integrand_G_R;
        integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_r, &neval);

        integrand_g.function = &Integrand_G_I;
        integration_romberg(&integrand_g, z1, z2, 20, 1e-10, 1e-7, &g_i, &neval);

        complex<double> value = complex<double>(g_r, g_i);

        return value;
    }

    complex<double> G1B(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        complex<double> value;

        double Rt2 = rho * rho + a * a + (z - z0) * (z - z0);
        double Rt = sqrt(Rt2);
        double Rt5 = Rt * Rt2 * Rt2;
        complex<double> emjkRt = complex<double>(cos(k * Rt), -sin(k * Rt));

        value = emjkRt / Rt + 0.25 * a * a * rho * rho
                * complex<double>(3.0 - k * k * Rt2, 3.0 * k * Rt) * emjkRt / Rt5;

        return value;
    }

    complex<double> G2B(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        complex<double> value;
        double Rt2 = rho * rho + a * a + (z - z0) * (z - z0);
        double Rt = sqrt(Rt2);
        double Rt3 = Rt * Rt2;
        double Rt5 = Rt2 * Rt3;
        complex<double> emjkRt = complex<double>(cos(k * Rt), -sin(k * Rt));

        if(rho > a) {
            value = 0.125 * a * a * rho * (2.0 + a * a / rho / rho)
                    * complex<double>(3.0 - k * k * Rt2, 3.0 * k * Rt) * emjkRt / Rt5;
            value = value + 0.5 * a * a * complex<double>(1.0, k * Rt) * emjkRt / rho / Rt3;
            value = value + emjkRt / rho / Rt;
        }
        else if(rho < a) {
            value = -0.125 * rho * rho * rho
                    * complex<double>(3.0 - k * k * Rt2, 3.0 * k * Rt) * emjkRt / Rt5;
            value = value - 0.5 * rho * complex<double>(1.0, k * Rt) * emjkRt / Rt3;
        }
        else {
            value = 0.5 * emjkRt / a / Rt + 0.125 * a * a * a
                    * complex<double>(3.0 - k * k * Rt2, 3.0 * k * Rt) * emjkRt / Rt5;
        }

        return value;
    }

    complex<double> G3B(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        complex<double> value;
        double Rt2 = rho * rho + a * a + (z - z0) * (z - z0);
        double Rt = sqrt(Rt2);
        double Rt3 = Rt * Rt2;
        complex<double> emjkRt = complex<double>(cos(k * Rt), -sin(k * Rt));

        if(rho > a) {
            value = 0.125 * complex<double>(0.0, k * a * a * rho) * (2.0 + a * a / rho /rho)
                    * complex<double>(1.0, k * Rt) * emjkRt / Rt3;
            value = value + emjkRt / rho + 0.5 * complex<double>(0.0, k * a * a) * emjkRt / rho / Rt;
        }
        else if(rho < a) {
            value = -0.125 * complex<double>(0.0, k * rho * rho * rho)
                    * complex<double>(1.0, k * Rt) * emjkRt / Rt3;
            value = value - 0.5 * complex<double>(0.0, k * rho) * emjkRt / Rt;
        }
        else {
            value = 0.5 * emjkRt / a + 0.125 * complex<double>(0.0, k * a * a * a)
                    * complex<double>(1.0, k * Rt) * emjkRt / Rt3;
        }

        return value;
    }

    complex<double> G1b(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        double Rt = sqrt(rho * rho + a * a + (z - z0) * (z - z0));
        complex<double> value = complex<double>(cos(k * Rt), -sin(k * Rt)) / Rt;
        return value;
    }

    complex<double> G2b(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        double Rt = sqrt(rho * rho + a * a + (z - z0) * (z - z0));
        complex<double> value = complex<double>(cos(k * Rt), -sin(k * Rt)) / rho / Rt;
        return value;
    }

    complex<double> G3b(double rho, double a, double z, double z0)
    {
        double k = M_2PI;
        double Rt = sqrt(rho * rho + a * a + (z - z0) * (z - z0));
        complex<double> value = complex<double>(cos(k * Rt), -sin(k * Rt)) / rho;
        return value;
    }

    complex<double> G4b(double rho, double a, double z, double z1, double z2)
    {
        double gt_r, gt_i;
        int neval;

        double k = M_2PI;

        INTEGRAND integrand_gt;

        Integrand_Parameter p;
        p.rho = rho;
        p.a = a;
        p.k = k;
        p.z = z;

        integrand_gt.parameter = &p;

        integrand_gt.function = &Integrand_Gt_R;
        integration_romberg(&integrand_gt, z1, z2, 20, 1e-10, 1e-7, &gt_r, &neval);

        integrand_gt.function = &Integrand_Gt_I;
        integration_romberg(&integrand_gt, z1, z2, 20, 1e-10, 1e-7, &gt_i, &neval);

        complex<double> value = complex<double>(gt_r, gt_i);

        return value;
    }

    static double Integrand_G_R(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double R = sqrt(p.rho * p.rho + (p.z - x) * (p.z - x));
        double value = cos(p.k * R) / R;
        return value;
    }

    static double Integrand_G_I(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double R = sqrt(p.rho * p.rho + (p.z - x) * (p.z - x));
        double value = -sin(p.k * R) / R;
        return value;
    }

    static double Integrand_G1_R(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double R = sqrt(p.rho * p.rho + p.a * p.a + p.z * p.z - 2.0 * p.a * p.rho * cos(x));
        double value = cos(p.k * R) / R;
        return value;
    }

    static double Integrand_G1_I(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double R = sqrt(p.rho * p.rho + p.a * p.a + p.z * p.z - 2.0 * p.a * p.rho * cos(x));
        double value = -sin(p.k * R) / R;
        return value;
    }

    static double Integrand_G2_R(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double R = sqrt(p.rho * p.rho + p.a * p.a + p.z * p.z - 2.0 * p.a * p.rho * cos(x));
        double rho = sqrt(p.rho * p.rho + p.a * p.a - 2.0 * p.a * p.rho * cos(x));
        double cosbeta = (p.rho - p.a * cos(x)) / rho;
        double value = cosbeta * cos(p.k * R) / rho / R;
        return value;
    }

    static double Integrand_G2_I(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double R = sqrt(p.rho * p.rho + p.a * p.a + p.z * p.z - 2.0 * p.a * p.rho * cos(x));
        double rho = sqrt(p.rho * p.rho + p.a * p.a - 2.0 * p.a * p.rho * cos(x));
        double cosbeta = (p.rho - p.a * cos(x)) / rho;
        double value = -cosbeta * sin(p.k * R) / rho / R;
        return value;
    }

    static double Integrand_G3_R(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double R = sqrt(p.rho * p.rho + p.a * p.a + p.z * p.z - 2.0 * p.a * p.rho * cos(x));
        double rho = sqrt(p.rho * p.rho + p.a * p.a - 2.0 * p.a * p.rho * cos(x));
        double cosbeta = (p.rho - p.a * cos(x)) / rho;
        double value = cosbeta * cos(p.k * R) / rho;
        return value;
    }

    static double Integrand_G3_I(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double R = sqrt(p.rho * p.rho + p.a * p.a + p.z * p.z - 2.0 * p.a * p.rho * cos(x));
        double rho = sqrt(p.rho * p.rho + p.a * p.a - 2.0 * p.a * p.rho * cos(x));
        double cosbeta = (p.rho - p.a * cos(x)) / rho;
        double value = -cosbeta * sin(p.k * R) / rho;
        return value;
    }

    static double Integrand_G4_R(double x, void * params)
    {
        double value;
        int neval;

        Integrand_Parameter p = * ((Integrand_Parameter *) params);

        INTEGRAND integrand_g;

        Integrand_Parameter pp;
        pp.rho = p.rho;
        pp.a = p.a;
        pp.k = p.k;
        pp.z = p.z - x;

        integrand_g.parameter = &pp;

        integrand_g.function = &Integrand_G1_R;
        integration_romberg(&integrand_g, 0.0, M_2PI, 20, 1e-10, 1e-7, &value, &neval);

        return value;
    }

    static double Integrand_G4_I(double x, void * params)
    {
        double value;
        int neval;

        Integrand_Parameter p = * ((Integrand_Parameter *) params);

        INTEGRAND integrand_g;

        Integrand_Parameter pp;
        pp.rho = p.rho;
        pp.a = p.a;
        pp.k = p.k;
        pp.z = p.z - x;

        integrand_g.parameter = &pp;

        integrand_g.function = &Integrand_G1_I;
        integration_romberg(&integrand_g, 0.0, M_2PI, 20, 1e-10, 1e-7, &value, &neval);

        return value;
    }

    static double Integrand_Gt_R(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double Rt = sqrt(p.rho * p.rho + p.a * p.a + (p.z - x) * (p.z - x));
        double value = cos(p.k * Rt) / Rt;
        return value;
    }

    static double Integrand_Gt_I(double x, void * params)
    {
        Integrand_Parameter p = * ((Integrand_Parameter *) params);
        double Rt = sqrt(p.rho * p.rho + p.a * p.a + (p.z - x) * (p.z - x));
        double value = -sin(p.k * Rt) / Rt;
        return value;
    }
}
